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Quantum  tunnelling  is  a  phenomenon  in  which  a  quantum  state  traverses  energy  barriers 
higher  than  the  energy  of  the  state  itself.  Quantum  tunnelling  has  been  hypothesized  as  an 
advantageous  physical  resource  for  optimization  in  quantum  annealing.  However, 
computational  multiqubit  tunnelling  has  not  yet  been  observed,  and  a  theory  of  co-tunnelling 
under  high-  and  low-frequency  noises  is  lacking.  Here  we  show  that  8-qubit  tunnelling  plays 
a  computational  role  in  a  currently  available  programmable  quantum  annealer.  We  devise  a 
probe  for  tunnelling,  a  computational  primitive  where  classical  paths  are  trapped  in  a  false 
minimum.  In  support  of  the  design  of  quantum  annealers  we  develop  a  nonperturbative 
theory  of  open  quantum  dynamics  under  realistic  noise  characteristics.  This  theory  accurately 
predicts  the  rate  of  many-body  dissipative  quantum  tunnelling  subject  to  the  polaron 
effect.  Furthermore,  we  experimentally  demonstrate  that  quantum  tunnelling  outperforms 
thermal  hopping  along  classical  paths  for  problems  with  up  to  200  qubits  containing  the 
computational  primitive. 
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^riV  uantum  annealing1-5  is  a  technique  inspired  by  classical 
I  1  simulated  annealing6  that  aims  to  take  advantage  of 
quantum  tunnelling.  In  classical  cooling  optimization 
algorithms  such  as  simulated  annealing,  the  initial  temperature 
must  be  high  to  overcome  tall  energy  barriers.  As  the  algorithm 
progresses,  the  temperature  is  gradually  lowered  to  distinguish 
between  local  minima  with  small  energy  differences.  This  causes 
the  stochastic  process  to  freeze  once  the  thermal  energy  is  lower 
than  the  height  of  the  barriers  surrounding  the  state.  In  contrast, 
quantum  tunnelling  transitions  are  still  present  even  at  zero 
temperature.  Therefore,  for  some  energy  landscapes,  one  might 
expect  that  quantum  dynamical  evolutions  can  converge  to  the 
global  minimum  faster  than  the  corresponding  classical  cooling 
process. 

The  goal  of  quantum  annealing  is  to  find  low-energy  states  of  a 
‘problem  Hamiltonian’ 

tfp  =  -  h w 

fi  fiV 

where  the  Pauli  matrices  correspond  to  spin  variables 
with  values  { ±  1}.  The  local  fields  { h ^  and  couplings  {J^v} 
define  the  problem  instance.  Quantum  annealing  is  characterized 
by  evolution  under  the  Hamiltonian 

H0(s)=A(s)HD+£(s)HP,  (2) 

where  HD  =  —  ^  cr*.  The  annealing  parameter  s  slowly 
increases  from  0  to  1  throughout  the  annealing  time  tqa. 
Initially,  A(0)  B( 0).  With  increasing  s,  A(s)  monotonically 

decreases  to  0  for  s=  1,  whereas  B(s)  increases. 

The  performance  of  chips  designed  to  implement  quantum 
annealing  using  superconducting  electronics  has  been  studied  in  a 
number  of  recent  lines  of  work7-24.  Given  a  quantum  annealer 
operating  at  finite  temperature,  noise  and  dissipation  strengths, 
does  it  utilize  tunnelling  or  thermal  activation  for  computation? 
Here  we  address  precisely  this  question.  We  introduce  a  16 -qubit 
probe  for  tunnelling,  a  computational  primitive  where  classical 
paths  are  trapped  in  a  false  minimum.  We  present  a 
nonperturbative  theory  of  multiqubit  tunnelling,  which  takes  into 
account  both  high-  and  low-frequency  noises.  To  distinguish 
between  tunnelling  and  thermal  activation,  we  study  the  thermal 
dependence  of  the  probability  of  success  for  the  computational 
primitive.  Thermal  activation  shows  an  increasing  probability  of 
success  with  increasing  temperature,  as  expected.  Multiqubit 
tunnelling,  on  the  other  hand,  shows  a  decreasing  probability  of 
success  with  increasing  temperature,  both  in  theory  and 
experiment.  Finally,  we  study  a  generalization  of  the 
computational  primitive  to  a  larger  number  of  qubits  that  contain 
the  same  ‘motif  multiple  times.  Quantum  tunnelling  outperforms 
thermal  hoping  for  these  problems  under  similar  parameters. 


Figure  1  |  Graph  of  the  tunnelling  probe  Hamiltonian.  The  16  qubits  are 
coupled  ferromagnetically  with  J  =  1  (lines).  The  applied  fields  are 
0<hL<J/2  (hR=  -1)  for  the  left  (right)  qubit  cell.  The  symmetry  and 
strong  intracell  ferromagnetic  coupling  makes  each  8-qubit  cluster  evolve 
together. 


Results 

Quantum  tunnelling  probe.  We  now  describe  a  probe  for  com¬ 
putational  tunnelling:  a  non- convex  optimization  problem 
consisting  of  just  one  global  minimum  and  one  false  (local) 
minimum.  For  concreteness,  we  use  the  problem  Hamiltonian 
depicted  in  Fig.  1,  implementable  in  a  D-Wave  Two  quantum 
annealer  (a  description  of  the  flux  qubits  used  in  D-Wave  Two  is 
given  in  the  Supplementary  Note  1).  The  problem  Hamiltonian 
consists  of  two  cells,  left  and  right,  each  with  n  —  8  qubits.  The 
local  fields  0</zL<0.5  and  hR  —  —  1  are  equal  for  all  the  spins 
within  each  cell,  and  all  the  couplings  J—\  are  ferromagnetic.  The 
spins  within  each  cell  tend  to  move  together  as  clusters  because  of 
symmetry  and  the  strong  intracell  ferromagnetic  coupling  energy. 
We  choose  |/zR|  >  \hL\  so  that  in  the  low-energy  states  of  HP  the 
right  cluster  is  pointing  along  its  own  local  field  as  seen  in  Fig.  1. 
The  difference  in  energy  of  the  states  with  opposite  polarization  in 
the  left  cluster  is  n(J—2hL).  Choosing  hL<J/ 2  =  0.5,  the  global 
minimum  corresponds  to  both  clusters  having  the  same  orienta¬ 
tion,  while  in  the  false  minimum  they  have  opposite  orientations. 

We  can  gain  an  intuitive  understanding  of  the  effective  energy 
landscape  if  we  represent  each  qubit  by  a  mean  field  spin  vector 
in  the  xz  plane.  Denote  by  6 ^  the  angle  of  the  spin  vector  for  qubit 
/i  with  the  x  quantization  axis.  We  assume  that  all  the  qubits  in 
the  left  (right)  cluster  have  the  same  angle  6L  ( 0R ).  This 
assumption  is  based  on  symmetry  and  the  strong  intracluster 
ferromagnetic  energy.  The  resulting  energy  potential  can  also  be 
derived  using  more  formal  methods,  such  as  the  Villain 
representation  (see  Supplementary  Note  2  and  ref.  25).  Figure  2 
plots  the  effective  energy  potential  for  the  left  cluster  as  a  function 
of  0L  with  hL  ==  0.44.  At  the  beginning  of  the  annealing  process 
£(s)/A(s)  C  1,  and  we  have  (crp  ~  h^B(s)/A(s)  (the  coupling 
terms  are  quadratic  in  the  z  polarizations  (crp  and  therefore 
negligible  at  this  point).  As  hL  and  hR  have  opposite  signs,  so  will 
the  z-projections  of  spins  in  the  two  clusters  early  in  the 
evolution.  To  escape  this  path  classically  all  spins  in  the  left 
cluster  must  flip  sign,  which  requires  traversing  an  energy  barrier. 
The  barrier  peak  corresponds  to  zero  total  z-polarization  of  the 
left  cluster.  Therefore,  the  barrier  grows  with  the  ferromagnetic 
energy  of  the  cluster  (n/2)2/.  The  barrier  height  is  much  greater 
than  the  residual  energy,  which  grows  with  n(J  —  2 hL). 


0.10 


Figure  2  |  Effective  energy  potential  using  fiL  =  0.44.  The  mean  field 
potential  is  plotted  versus  annealing  parameter  s  and  tilt  angle  0L  of  each 
spin  vector  in  the  left  cluster.  The  red  line  corresponds  to  a  path  that  starts 
in  the  initial  global  minimum  and  follows  the  instantaneous  local  energy 
minimum.  A  second  local  minimum  (dashed  blue  line)  forms  at  the 
bifurcation  point  s  =  0.18.  The  global  minimum  is  found  in  this  second  path 
after  s  =  0.24  (dashed  to  continuous  blue  line).  To  reach  this  global 
minimum  the  system  state  has  to  traverse  the  energy  barrier  between  them 
(dashed  green  line),  either  by  thermal  activation  or  by  quantum  tunnelling. 
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Quantum  mechanically,  if  hL<J/ 2  the  system  evolution  goes 
through  an  ‘avoided-crossing’  where  the  two  lowest  eigenstates 
Ei(s)  and  E0(s)  approach  closely  to,  and  then  repel  from,  each 
other  (see  inset  in  Fig.  3).  Higher-energy  states  remain  well 
separated  during  the  evolution.  This  level  repulsion  occurs 
because  of  the  collective  tunnelling  of  qubits  in  the  left  cluster 
between  the  opposite  z  polarizations.  At  the  point  where  the  gap 
hQi0(s)  =  Ei(s)  —  E0(s)  reaches  its  minimum,  the  corresponding 
adiabatic  eigenstates  are  formed  by  the  symmetric  and  antisym¬ 
metric  superpositions  of  the  cluster  orientations.  The  size  of  the 
minimum  gap  varies  with  hL  as  seen  in  Fig.  3.  The  position  of  the 
avoided-crossing  can  be  estimated  to  occur  at  the  point  where 
1(^)1  ~  2 hL/J  and  moves  towards  s=  1  as  approaches  J/2. 
Note  that  for  hL  =  J/2  the  residual  energy  n(J  —  2 hL)  vanishes  and 
the  final  ground  space  is  degenerate.  There  is  no  avoided-crossing 
when  hL>J/ 2  (Supplementary  Fig.  1). 

Characterization  of  noise  and  dissipation.  Under  realistic  con¬ 
ditions,  a  quantum  annealer  can  be  strongly  influenced  by  coupling 
to  the  environment.  We  introduce  a  detailed  phenomenological 
open  quantum  system  model  based  on  single-qubit  measurable 
noise  parameters.  We  shall  assume  that  each  qubit  is  coupled  to  its 
own  environment  with  an  independent  noise  source.  In  the  con¬ 
crete  case  studied  here,  this  is  consistent  with  experimental  data9, 
and  the  coupling  of  the  environment  to  each  flux  qubit  is 
proportional  to  a  oz  qubit  operator  (flux  fluctuations). 

In  the  analysis  of  the  transitions  between  the  states  we  start 
from  the  initial  (gapped)  stage  when  the  instantaneous  energy  gap 
hQi0(s)  between  the  two  lowest  eigenstates  |i /r0(s)),  |*Ai (s))  is 
sufficiently  large  compared  with  the  linewidth  hW.  Then,  the 
coupling  to  the  environment  can  be  treated  as  a  perturbation  and 
the  transition  rate  between  these  states  is  given  by  Fermi’s  golden 
rule  ri_>0(s)  ~  a(s)S(£2io(s))/fi2,  where  S(co)  is  the  noise 
spectral  density  (see  Methods).  Here 

2n  I  |2 

a(s)  =  J2\(Us)K\^1(s))\  (3) 

0=1 

is  a  sum  of  (squared)  transition  matrix  elements  between  the  two 
eigenstates. 


Figure  3  |  Quantum  energy  gap.  Inset  shows  the  gap  M210  =  F|(s)  -  E0(s) 
versus  s,  using  /il  =  0.44.  The  dashed  line  is  the  gap  in  the  diabatic 
(pointer)  basis.  In  the  main  plot,  the  minimum  gap  decreases  with  hL.  The 
horizontal  green  dashed  line  (324  MHz)  corresponds  to  15.5  mK,  the  lowest 
temperature  in  our  experiments.  The  lower  inset  shows  the  spin 
configurations  of  the  two  lowest  eigenstates  at  the  end  of  the  annealing. 


In  the  minimum  gap  region,  the  (squared)  matrix  element  a(s) 
for  the  transition  rate  is  large,  and  the  system  is  thermalized 
(Fig.  4).  More  precisely,  we  have  1  /tqa,  where  the 

inverse  of  the  annealing  time  1  /tqa  is  an  approximation  for 
the  annealing  rate.  The  ground-state  population  is  given  by  the 
Boltzmann  distribution  at  the  experimental  temperature. 

After  the  avoided-crossing  region  (at  s  =  0.255)  we  observe  a 
steep  exponential  fall- off  of  the  matrix  element  a(s)  with  s, 
eventually  causing  multiqubit  freezing  (Fig.  4).  Multiqubit 
freezing  is  quite  distinct  from  single-qubit  freezing.  Single-qubit 
tunnelling10  decays  slowly  as  the  magnitude  of  the  transverse  field 
A(s)  decreases.  The  multiqubit  transition  rate,  however,  decays 
exponentially  fast  (see  inset  of  Fig.  4).  This  is  due  to  the 
increasing  effective  barrier  width  (Fig.  2),  which  results  in  an 
exponential  decrease  in  quantum  tunnelling  and  a  slowdown  of 
the  transition  rate  Formally,  the  barrier  width  corresponds 

to  the  Hamming  distance 

2  n  2 

fcW  =  E|^o|^o)-(^|^|^)|  /4  (4) 

0=1 

between  the  opposite  z  orientations  of  the  left  cluster  in  the  two 
lowest-energy  eigenstates.  The  exponential  sensitivity  of 
multiqubit  tunnelling  to  the  width  or  Hamming  distance  h(s)  is 
the  cause  of  the  exponential  decay  of  the  matrix  element  a(s),  and 
of  the  multiqubit  freezing. 

We  distinguish  a  slowdown  phase  (roughly,  0.1  <  tqaTi^0<  10) 
and  a  frozen  phase  (^ar1^o<0-l).  In  the  frozen  phase,  there  are 
no  dynamics.  Part  of  the  system  population  remains  trapped  in 
the  excited  state  |i/q(s))  corresponding  to  the  false  minimum  of 
the  effective  potential  until  the  end  of  the  quantum  annealing 
process  (Fig.  4). 


Figure  4  |  Multiqubit  freezing.  Solid  lines  correspond  to  the  modelled 
population  pm  of  the  lowest-energy  eigenstate  along  the  quantum  annealing 
process  using  hL  =  0.44  at  15.5  mK  (top  line)  and  35  mK  (bottom  line). 
Dashed  lines  correspond  to  the  thermal  equilibrium  population  pB.  In  the 
thermalization  phase  (red)  the  transition  rate  is  fast  and  the  population 
remains  close  to  thermal  equilibrium.  As  the  multiqubit  energy  barrier 
increases,  the  transition  rates  are  exponentially  reduced  with  s,  as  shown  in 
the  inset.  We  define  the  slowdown  regime  (blue)  as  n ^ 0 <  1 0  and  the 

frozen  regime  (grey)  as  tqar1_>0<0.1.  Comparing  the  data  at  15.5  and 
35  mK,  we  see  a  small  change  in  the  transition  rate  relative  to  the  larger 
change  in  the  thermal  equilibrium  ground-state  population.  Therefore,  the 
probability  of  success  is  lower  at  higher  temperature. 
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When  the  energy  gap  is  similar  to  (or  smaller  than)  the  noise 
linewidth  W,  the  environment  cannot  be  treated  as  a  perturba¬ 
tion.  We  develop  a  multiqubit  nonperturbative  analysis  in  the 
spirit  of  the  Non-interacting  Blip  Approximation  (NIBA)26  that 
covers  all  quantum  annealing  (QA)  stages.  In  the  slowdown 
phase,  when  the  Hamming  distance  approaches  its  maximum 
value,  h  ~  n,  the  instantaneous  decay  rate  of  the  first  excited  state 
takes  the  form 


Ti 


40  =  J  dr  e‘^-h(kPT+(Wz)2/2)  |y£csc/,  (T  nL) 


P/n  _ 


D(  t),  (5) 


where  W  and  ep  are  the  linewidth  and  reconfiguration  energy  from 
the  low-frequency  noise,  r\  and  1  /tc  are  the  high-frequency  Ohmic 
noise  coupling  and  cutoff  and  /?  is  the  inverse  temperature  h/kBT. 
The  dependence  on  the  annealing  parameter  s  is  implicit.  The 
factor  D(t)  is  related  to  the  tunnelling  permeability  of  the  potential 
barrier  in  Fig.  2  (similar  to  the  coefficient  a).  It  has  the  expression 

D(t,  s)  =  [(ep(s)c+  (s)  -  ic _  (s)<9t£(t,  s))2  +  a(s)dTTf  (t,  s)j  e~tc^d^ 

(6) 

where 


Zyy  M  _ 


'oo  =  (yoo|^|</vOo),  v,y'  =  0,1 

1  2n  /  X 

C±(5)  (5)  ±Z«°(s)) 

p=\ 

<>M=jE  (z”w)  -(z»”w) 


M=1  L 


C(t,5)  =  fep(5)l  + 


(W(5)t) 


(7) 


Equation  (5)  describes  collective  tunnelling  of  the  left  qubit 
cluster  assisted  by  the  environment.  The  crucial  difference  from 
the  single- qubit  theory27,28  is  that  the  parameters  of  the 
environment  in  the  transition  rate  are  rescaled  by  the  barrier 
width  or  Hamming  distance  h(s).  The  effective  low-frequency 
noise  linewidth  is  /z1/2(s)  W(s),  the  reconfiguration  energy  is 
h(s)ep(s)  and  the  Ohmic  coefficient  is  h(s)rj(s).  This  is  important 
at  the  late  stages  of  quantum  annealing  when  h  ~  n  1. 


Comparison  of  NIBA  with  data.  We  observe  a  very  close  cor¬ 
respondence  between  the  results  of  the  analysis  with  the  NIBA 
Quantum  Master  Equation  for  the  dressed  cluster  states  and  the 
D-Wave  Two  data  displayed  in  Fig.  5,  which  shows  the  depen¬ 
dence  of  the  ground-state  population  on  hL.  We  emphasize  that 
for  NIBA  (and  the  standard  Redfield  equation  with  Ohmic 
50h(ft))>  see  Methods)  we  do  not  have  any  parameter  fitting:  the 
parameters  are  obtained  from  experiments  (see  Methods).  The 
success  probability  of  quantum  annealing  is  (roughly)  determined 
by  the  thermal  equilibrium  ground-state  population  during  the 
slowdown  phase  (Fig.  4).  When  the  temperature  grows,  the 
ground-state  population  decreases  appreciably,  while  the  transi¬ 
tion  rate  changes  little.  Consequently,  quantum  mechanically,  the 
probability  of  success  decreases  with  increasing  temperature,  as 
seen  in  Fig.  6,  for  sufficiently  big  gaps.  Figure  6  shows  the 
dependence  of  the  ground-state  population  with  temperature. 

For  hL  closer  to  the  degeneracy  value  hL  =  //2,  the  minimum 
gap  becomes  smaller,  as  seen  in  Fig.  3.  Where  Q i0  <C  W, 
the  adiabatic  basis  of  the  instantaneous  multiqubit  states 
{ |i/^0(^))  ?  |i/q(s))}  loses  its  physical  significance.  Because  the 
coupling  to  the  bath  is  relatively  strong  here,  the  system  quickly 
approaches  the  states  corresponding  to  predominantly  opposite 
cluster  orientations,  similar  to  diabatic  states  (see  inset  of  Fig.  3). 


Figure  5  |  Probability  of  success  versus  fiL.  We  plot  the  probability  of 
measuring  the  final  ground  state  for  different  values  of  hL.  The  physical 
temperature  for  D-Wave  is  15.(5)  mK,  and  the  annealing  time  is  20  |is.  Both 
Redfield  and  NIBA  use  only  measured  parameters  (no  fitting).  SVMC  uses  an 
algorithmic  temperature  equal  to  the  physical  temperature  and  128,000 
sweeps,  as  explained  in  the  text.  Error  bars  (s.e.m.)  are  smaller  than  markers. 


Transitions  between  these  states,  also  called  pointer  states29, 
occur  at  a  much  slower  rate  as  a  consequence  of  the  polaronic 
effect.  As  a  result,  for  sufficiently  small  mininum  gaps  the 
multiqubit  freezing  starts  before  the  avoided- crossing  and  the 
success  probability  increases  with  temperature12.  This  is  captured 
by  the  multiqubit  NIBA  equation,  but  not  by  the  standard 
Redfield  Quantum  Master  Equation. 

Spin  vector  Monte  Carlo.  We  want  to  distinguish  quantum 
tunnelling  from  thermal  activation  along  classical  paths  of 
product  states  (which  preclude  multiqubit  tunnelling).  To  give  a 
more  precise  description  of  the  classical  paths  of  product  states, 
let  each  qubit  be  represented  by  a  mean  field  spin  vector  in 
the  xz  plane  and  denote  by  6 ^  the  angle  of  the  spin  vector  for 
qubit  with  the  x  quantization  axis,  as  before.  A  classical  path 
(red  line  in  Fig.  2)  that  follows  the  local  minimum  of  the  effective 
energy  potential  gets  trapped  in  a  false  minimum  and  fails  to 
solve  the  corresponding  optimization  problem,  as  explained 
above.  In  the  absence  of  quantum  tunnelling,  the  global  mini¬ 
mum  could  be  reached  through  thermal  excitations  for  over-the- 
barrier  escape  from  the  false  minimum.  This  thermal  activation 
results  in  an  increasing  probability  of  success  with  rising 
temperature. 

This  intuition  is  supported  by  spin  vector  Monte  Carlo 
(SVMC),  a  numerical  algorithm  consisting  of  thermal  Metropolis 
updates  of  the  spin  vectors.  SVMC  was  introduced  recently  in 
ref.  20  and  studied  in  related  lines  of  work21,23.  The  dynamics  are 
constrained  to  spin-vector  product  states,  with  one  spin  vector 
per  qubit.  For  a  given  Hamiltonian  H0(s),  we  denote  the 
corresponding  energy  by  Es(6i>...  ,0nq),  where  nq  is  the  number 
of  qubits.  The  evolution  consists  of  a  sequence  of  sweeps  along 
the  Hamiltonian  path  {H0(s)}.  In  each  sweep,  a  Monte  Carlo 
update  is  proposed  for  each  qubit  in  two  steps.  First,  a  new  angle 
6'^ is  drawn  from  the  uniform  distribution  in  [0,  2n].  Second, 
the  spin  vector  for  qubit  fi  is  updated  ^according  to  the 

Metropolis  rule  for  the  energy  difference 

D  =  £s(...  ,0/  ...) -£,(•••  A,  •••)• 
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Figure  6  |  Probability  of  success  versus  temperature  for  hL  =  0.44.  The 

decreasing  probability  with  increasing  temperature  is  only  matched  with 
theories  based  on  quantum  tunnelling.  This  is  the  opposite  tendency  to 
thermal  activation  (SVMC).  The  annealing  time  is  20  |is.  Both  Redfield  and 
NIBA  use  only  measured  parameters  (no  fitting).  SVMC  uses  an  algorithmic 
temperature  equal  to  the  physical  temperature  and  128,000  sweeps.  In  this 
temperature  range  the  lowest  two  states  (the  double-well  potential) 
account  for  all  the  probability  (0.9998  for  D-Wave,  0.99998  for  SVMC). 
Error  bars  (s.e.m.)  are  smaller  than  markers. 

That  is,  the  move  is  always  accepted  if  D  is  negative,  and  with 
probability  given  by  the  Boltzmann  factor  exp(  —  D/kBT)  if  D  is 
positive. 

The  initial  state  is  chosen  to  be  the  global  minimum  of  the 
transverse  field.  For  low  T  and  sufficient  sweeps,  the  evolution 
proceeds  along  the  false  minima  path  of  Fig.  2.  This  numerical 
method  allows  us  to  study  thermal  hopping  between  the  classical 
paths.  To  check  this  correspondence,  we  studied  the  height  of  the 
energy  barrier  obtained  from  Kramers’  theory  applied  to  SVMC. 
For  the  effective  potential  at  a  fixed  value  of  s,  we  initialized  the 
spin  vector  state  at  a  local  minima,  and  watch  for  Kramers  events. 
A  Kramers  event  corresponds  to  the  arrival  at  the  other  minima 
under  thermal  activation.  According  to  Kramers’  theory,  the 
dependence  on  temperature  for  the  expected  number  of  sweeps 
necessary  for  a  Kramers  event  follows  the  formula  exp  (A  17/ T), 
where  AC  is  the  height  of  the  energy  barrier.  We  extract  the 
energy  barrier  by  fitting  the  curve  of  the  average  number  of 
sweeps  for  different  T.  We  find  that  this  matches  almost  exactly 
the  energy  barrier  height  in  Fig.  2  for  different  values  of  s  (Fig.  7). 

A  disadvantage  of  SVMC  as  outlined  above  and  introduced  in 
ref.  20  is  that  there  is  no  natural  choice  to  relate  the  number  of 
sweeps  to  the  physical  evolution  time.  As  in  other  lines  of  work, 
we  will  choose  the  number  of  sweeps  to  obtain  a  good  correlation 
with  the  probability  of  success  of  the  D-Wave  chip  for  a 
benchmark  of  random  Ising  models  with  binary  couplings 
/^vg{1,  —1}  (refs  15,18,20,23).  This  will  allow  us  to 

phenomenologically  correlate  the  number  of  sweeps  to  physical 
time.  We  set  the  algorithmic  temperature  of  SVMC  to  be  the 
same  as  the  physical  temperature  because  we  are  interested  in  the 
dependence  of  the  success  probability  with  temperature.  There 
are  no  important  differences  for  the  correlation  with  other 
temperature  choices.  The  correlation  with  the  random  Ising 
benchmark  for  128,000  sweeps  (Fig.  8)  is  0.92,  and  the  residual 
probabilities  Psvmc  —  ^D-wave  have  a  mean  of  0.05  and  a  s.d.  of 
0.12.  This  is  consistent  with  the  best  values  found  over  a  wide 
range  of  parameters.  We  therefore  use  128,000  sweeps  at  15  mK 
as  our  reference  parameters  for  SVMC. 


Figure  7  |  Analysis  of  the  activation  energy  for  Kramer's  scape  for 

SVMC.  The  main  figure  shows,  in  a  semilog  scale,  the  average  number  of 
sweeps  as  a  function  of  temperature.  We  plot  lines  for  different  points  in 
the  annealing  schedule,  from  s  =  0.217  (dark  blue,  see  legend)  to  s  =  0.265 
(green).  Error  bars  correspond  to  the  s.e.m.  The  embedded  figure  shows 
the  activation  energy  (red  dots)  and  the  semiclassical  energy  barrier  (blue). 
There  is  a  good  correspondence  between  SVMC  and  the  effective  energy 
potential. 

Figure  6  confirms  the  thermal  activation  in  SVMC.  This  is 
opposite  to  both  open  quantum  system  theory  and  experiments 
with  the  D-Wave  chip,  which  show  a  reduction  in  the  probability 
of  success  with  rising  temperature,  as  explained  above. 
Furthermore,  Fig.  5  shows  that  the  probability  of  success  for 
SVMC  is  lower  than  the  probability  of  success  for  D-Wave  and 
open  system  quantum  models. 

Larger  problems.  A  generalization  of  the  16-qubit  problem  to  a 
larger  number  of  qubits  is  achieved  by  studying  problems  that 
contain  the  same  ‘motif  (Fig.  1)  multiple  times  within  the 
connectivity  graph  (Fig.  9).  The  success  probabilities  for  up  to  200 
qubits  are  shown  in  Fig.  10.  We  fit  the  average  success  probability 
as  P(nq)ocexp(  —  anq),  where  nq  is  the  number  of  qubits. 
The  fitting  exponent  a  for  the  D-Wave  Two  data  is 
(1.1  ±  0.05)  x  10  _2,  while  the  fitting  exponent  for  the  SVMC 
numerics  is  (2.8  ±0.17)  x  10  ~2.  We  conclude  that,  for  instance 
with  multiqubit  quantum  tunnelling,  the  D-Wave  Two  processor 
returns  the  solution  that  minimizes  the  energy  with  consistently 
higher  probability  than  physically  plausible  models  of  the 
hardware  that  only  employ  product  states  and  do  not  allow  for 
multiqubit  tunnelling  transitions. 

Discussion 

The  role  of  multiqubit  tunnelling  as  a  computational  resource  is  an 
open  problem  of  active  research.  Nevertheless,  it  is  instructive  to 
consider  some  plausible  estimates  for  the  case  of  minimization  of  an 
Ising  problem  that  contains  pairwise  interactions  between  all  qubits. 
A  way  to  think  of  multiqubit  tunnelling  as  a  computational  resource 
is  to  regard  it  as  a  form  of  large  neighbourhood  search.  Collective 
tunnelling  transitions  involving  K  qubits  explore  a  K  variable 
neighbourhood,  and  there  is  a  combinatorial  number  of  such 
neighbourhoods.  Using  standard  resources,  the  cost  of  exhaustively 
searching  on  a  Hamming  ball  of  binary  strings  of  radius  K  is 
Y^f=  i  (”)>  which  is  bounded  from  below  by  (”),  where  n  is  the 
total  number  of  qubits.  Therefore,  for  K  <C  n,  the  cost  is 
~  nK /K\  exp (K log («)).  As  one  can  see,  the  exponent  in  K 

is  very  steep  (log  n).  On  the  other  hand,  for  K~n/2  the  exponent  is 
that  of  the  exhaustive  search  in  the  entire  n— bit  string  space  (2"). 
We  now  compare  this  with  the  tunnelling  rate  across  a  barrier  of 
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15  mK  16-K  sweeps  1 5  mK  32-K  sweeps  15  mK  64-K  sweeps  15  mK  128-K  sweeps  15  mK  256-K  sweeps  15  mK  51 2-K  sweeps  15  mK  1 ,024-K  sweeps 


10  mK  16-K  sweeps  10  mK  32-K  sweeps  1 0  mK  64-K  sweeps  10  mK  128-K  sweeps  10  mK  256-K  sweeps  10  mK  51 2-K  sweeps  10  mK  1 ,024-K  sweeps 


Figure  8  |  Scatter  plots  showing  the  correlation  of  D-Wave  Two  data  and  SVMC.  Correlations  for  the  random  Ising  benchmark  for  different  algorithmic 
temperatures  (in  mK)  and  number  of  sweeps.  We  will  use  the  parameters  7=  15  mK  and  sweeps  =  128,000  for  SVMC  in  the  rest  of  the  paper. 


Figure  10  |  Success  probability  for  a  glass  of  clusters  as  a  function  of  the 
number  of  qubits.  We  fit  the  mean  probability  of  success  P(nq)  oc  exp(  -  a/iq) 
as  a  function  of  the  number  of  qubits  nq  (dashed  lines).  The  fitting  exponent  a 
for  the  D-Wave  Two  data  is  (1.1  ±  0.05)  x  10  ~  2  while  the  fitting  exponent  for 
the  SVMC  numerics  is  (2.8  ±0.17)  xIO-2  The  error  estimates  for  the 
exponents  (s.e.m.)  are  obtained  by  bootstrapping. 


Figure  9  |  Larger  problems  that  contain  the  tunnelling  probe  'motif  as 
subproblems.  As  in  Fig.  1,  the  black  (grey)  cluster  has  a  strong  hR  =  - 1 
(weak  /\  =  0.44)  local  field.  The  black  clusters  are  connected  in  a  glassy 
manner  to  make  the  problem  less  regular:  all  connections  between  any  two 
neighbouring  black  clusters  are  set  randomly  to  either  —  1  or  +1.  The  —  1 
connections  are  depicted  in  blue. 


width  K.  It  is  given  by  exp(  —  cK)  for  some  constant  c.  If  K  is 
smaller  than  n  then,  as  n  increases,  exp(cTC)  <C  exp(TC  log  («))  and 
tunnelling  can  be  faster  than  large  neighbourhood  search.  On  the 
other  hand,  in  many  problems  similar  to  the  two-cluster  problem 
the  barrier  width  is  K  in  0(n).  In  these  cases  the  tunnelling  rate  can 
still  be  exp  (cri)  <C  2n  (refs  30-35).  Therefore,  again  tunnelling  can 
still  provide  a  dramatically  faster  search  option. 

We  find  that  the  current-generation  D-Wave  Two  annealer 
enables  tunnelling  transitions  involving  at  least  8  qubits.  It  will  be 
an  important  future  task  to  determine  the  maximal  K  attainable 
by  current  technology  and  how  large  it  can  be  made  in  next 
generations.  The  multiqubit  polaronic  quantum  master  equation 
presented  here  lays  down  the  theory  to  answer  this  question.  It 
guides  the  design  of  next-generation  architectures  and  helps  to 
understand  for  which  computational  problems  quantum- 
enhanced  optimization  may  offer  an  advantage.  The  larger  the 
K  the  easier  it  should  be  to  translate  the  quantum  resource 
‘K- qubit  tunnelling  into  a  possible  computational  speedup.  We 
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want  to  emphasize  that  this  paper  does  not  claim  to  have 
established  a  quantum  speedup.  To  this  end,  one  would  have  to 
demonstrate  that  no  known  classical  algorithm  finds  the  optimal 
solution  as  fast  as  the  quantum  process.  To  establish  such  an 
advantage  it  will  be  important  to  study  to  what  degree  collective 
tunnelling  can  be  emulated  in  classical  algorithms  such  as 
Quantum  Monte  Carlo  or  by  employing  cluster  update  methods. 
However,  the  collective  tunnelling  phenomena  demonstrated  here 
present  an  important  step  towards  what  we  would  like  to  call  a 
physical  speedup:  a  speedup  relative  to  a  hypothetical  version  of 
the  hardware  operated  under  the  laws  of  classical  physics. 


Methods 

Experimental  properties  of  the  noise.  The  properties  of  the  noise  are  determined 
by  the  noise  spectral  density  S(co),  which  is  characterized  by  single-qubit 
macroscopic  resonant  tunnelling  (MRT)  experiments  in  a  broad  range  of  biases 
(0.4  MHz-4  GHz)  and  temperatures  (21-38  mK)  for  tunnelling  amplitudes  of  a 
single  flux  qubit  below  1  MHz.  The  MRT  data  collected  are  surprisingly  well 
described28,36  by  a  phenomenological  ‘hybrid’  thermal  noise  model 
S(co)  =  Sif(co)  +  S0h(co).  Here  S0h(<w)  =  h2rj(De~(OXc / (l  —  e~h(°/kBT )  denotes  the 
high-frequency  part  and  has  Ohmic  form  with  dimensionless  coupling  rj  and  cutoff 
frequency  1/tc  (assumed  to  be  very  large).  The  low-frequency  part  Sif  is  of  the 
l//type37,  and  in  current  D-Wave  chips  this  noise  is  coupled  to  the  flux  qubit 
relatively  strongly.  Its  effect  can  be  described  with  only  two  parameters:  the  width 
W  and  the  Stokes  shift  ep  of  the  MRT  line27.  The  experimental  shift  value  is  related 
to  the  width  by  the  fluctuation-dissipation  theorem  (ep  =  h  W2/2kB  T)  and 
represents  the  reorganization  energy  of  the  environment.  The  values  of  the  noise 
parameters  measured  at  the  end  of  the  annealing  (s=  1)  for  the  D-Wave  Two  chip 
are  Wl(2n)  =  0.40(1)  GHz  and  r\  =  0.24(3). 
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